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The classical notion of a single-particle scalar distribution function or phase space density can 
be generalized to a matrix in order to accommodate superpositions of states of discrete quantum 
numbers, such as neutrino mass/flavor. Such a 'neutrino distribution matrix' is thus an appropriate 
construct to describe a neutrino gas that may vary in space as well as time and in which flavor 
mixing competes with collisions. The Liouville equations obeyed by relativistic neutrino distribution 
matrices, including the spatial derivative and vacuum flavor mixing terms, can be explicitly but 
elegantly derived in two new ways: from a covariant version of the familiar simple model of flavor 
mixing, and from the Klein-Gordon equations satisfied by a quantum 'density function' (mean value 
of paired quantum field operators). Associated with the latter derivation is a case study in how the 
joint position/momentum dependence of a classical gas (albeit with Fermi statistics) emerges from 
a formalism built on quantum fields. 
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I. INTRODUCTION 

The decoupling of neutrinos from dense nuclear matter 
occurs in a number of environments in which neutrino fla- 
vor mixing may also play an important role, including the 
early universe [HQ and core-collapse supernovae [!, 0, [j| . 
Treatment of the neutrinos' transition from diffusion to 
free streaming requires some sort of 'transport calcula- 
tion' embodying the principles of some kind of 'kinetic 
theory. ' 

Localized 'microscopic' quantum mechanical effects 
can be handled within the framework of classical kinetic 
theory [1, 0, In classical kinetic theory a single- 
particle distribution function /(i,x, p) quantifies the av- 
erage number dN of particles of a particular type, having 
spin degeneracy g and momenta within d 3 p of p, at po- 
sitions within d 3 x of x: 



dAT = /(i,x,p) 



gd 3 p 

(2tt) 3 



(1) 



The Boltzmann equation equates a collision integral C (/) 
to the rate of change df/dX of the average density of 
particles in classical phase space trajectories (cc(A), p(A)). 
Given the geodesic equations defining worldlines x(X), 



dx^ 

~d\ 

dp^ 

~d\ 



(2) 
(3) 



(here r M ^ p are the connection coefficients associated with 



the spacetime metric), the Boltzmann equation can be 
expressed 



Of 



df 



dp 1 



(4) 



On the left-hand side the Liouville operator acts upon 
/. The collision integral on the right-hand side repre- 
sents the phase space density (rate per space volume 
per momentum space volume) of isolated 'microscopic' 
or point-like scattering events between classical trajecto- 
ries. By reasonable extension, the collision integral also 
includes inherently quantum-mechanical processes affect- 
ing the population of classical phase space trajectories, 
such as particle decays, particle emission/absorption, and 
pair creation/annihilation. The restriction to isolated 
point-like transitions allows for insertion into C(/) of 
interaction rates computed (for instance) with the stan- 
dard methods of quantum field theory, together with fac- 
tors 1 ± / (with upper sign for bosons and lower sign 
for fermions) encoding the impact of quantum statistics 
upon available final-state phase space. 

However, neutrino flavor mixing is a 'macroscopic' 
quantum mechanical effect, requiring the evolution of 
amplitudes across time and/or distance scales compa- 
rable to scales characteristic of the total system under 
consideration; hence flavor mixing cannot in general be 
described by the scalar distribution function /(t,x, p) 
and the Boltzmann equation it obeys, these being con- 
cerned only with the evolution of particle number, and 
not quantum amplitudes and the evolution of phases. 

Classical transport being inadequate as a conceptual 
framework for handling flavor mixing, attention turns 
to a statistical treatment of a quantum gas. An ap- 
proach especially suited to cases of spatial homogene- 
ity begins with the definition of a quantum occupation 
number n{t, q), which quantifies the average number dN 
of particles of spin degeneracy g, occupying momentum 
eigenstates within d 3 q of q within the (effectively infi- 
nite) quantization volume V: 



dN = n(t, q) 



gVd 3 q 
(2tt) 3 ' 



(5) 
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The time derivative dn/dt is equal to an entity similar 
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to the collision integral in the Boltzmann equation, con- 
structed using transition rates between quantum states. 

Distinctions between quantum occupation numbers 
n(t, q) and classical distribution functions /(£,x, p) 
should be kept in mind. Because they specify the average 
populations of quantum states rather than classical tra- 
jectories, quantum occupation numbers are distinguished 
from classical distribution functions by an absence of 
spatial dependence, as required by the impossibility in 
quantum mechanics of simultaneous sharp specification 
of both position and momentum. Moreover, there is a 
subtle difference between the momenta p in f(t, x, p) and 
q in n(£,q). In the case of an occupation number, q is 
the eigenvalue of a momentum eigenstate. But in a clas- 
sical distribution function, p represents the momentum 
of a classical particle; considered as a limit of a quantum 
mechanical description, p might therefore be thought of 
as representing the centroid of a momentum-space wave 
packet — that is, p is the expectation value of a superpo- 
sition of momentum eigenstates q. 

Accommodation of flavor mixing in a neutrino gas 
was first contemplated [j| [l(| for the (homogeneous and 
isotropic) early universe, through the introduction of 
what might at first glance be thought of as a 'neutrino 
(quantum) occupation matrix' p(t, |q|). This is "a ma- 
trix in the space of neutrino species" [9| , whose diagonal 
elements are occupation numbers of the various neutrino 
species, and whose off-diagonal elements quantify the ex- 
tent to which neutrinos exist in superpositions of distinct 
species. Its introduction appears to have been motivated 
by the recognition that neutrino interaction amplitudes 
constitute a matrix with nonvanishing off-diagonal en- 
tries when written in terms of 'physical states' (neutrinos 
of definite mass) . The fact that the off-diagonal elements 
of a neutrino occupation matrix p(t, |q|) contain informa- 
tion on coherent superpositions is reminiscent of a density 
matrix, as is the fact that it obeys a Heisenberg-like equa- 
tion of motion (time derivative given by a commutator 
with a Hamiltonian) . Indeed it has sometimes been called 
a 'density matrix' in the literature (for instance, in Refs. 
H> EH)- However, it is perhaps better thought of as 
a 'matrix of densities' [l2| , since it is not a density ma- 
trix (as traditionally defined) for multiparticle neutrino 
states, in that its trace is not equal to unity. Hence its 
diagonal elements are n ot 'p robabilities' in the strictest 
(that is, absolute) sense [Ig[37j] [II]. 

In typical studies of the epoch of big-bang nucleosyn- 
thesis, however, the momenta apparently take on their 
classical significance — that is, there seems to be a (gener- 
ally tacit) assumption that the neutrino gas is described 
by a 'distribution matrix' pit, |p|) rather than an 'occu- 
pation matrix' p(t, |q|). This is because the expansion 
of the universe must be accounted for. Classical kinetic 
theory calculations in the context of the early universe — 
which do not involve flavor mixing — have long used the 
Boltzmann equation, or a momentum integral thereof. 
The redshift (or number dilution, in the momentum- 
integrated case) due to cosmological expansion results 



from nonvanishing connection coefficients in Eq. ([4]); 
see for instance Ref. [Hj]. The very reasonable, if un- 
remarked [j| [l(| , assumption seems to be that a replace- 
ment 

obtains, where the Hubble parameter H is the cosmo- 
logical expansion rate. That is, the total time derivative 
dp/dt of a quantum occupation matrix p(t, |q|) that sat- 
isfies a Heisenberg-like equation of motion somehow goes 
over to the action of the classical Liouville operator upon 
a distribution matrix p(t, |p|) that is classical in all but 
the discrete quantum numbers (e.g. flavor/mass) respon- 
sible for the matrix structure. The intuition behind this 
replacement is evidently similar to that which motivates 
the use of interaction rates computed with quantum field 
theory (plus quantum statistics) in the Boltzmann equa- 
tion's collision integral, as mentioned above in connection 
with classical kinetic theory. In particular, the 'Liouville 
replacement' of Eq. © and the calculation of interaction 
rates that go into collision integrals share the following 
feature: plane waves (momentum eigenstates), here la- 
beled by q, are taken as proxies for classical particles (or 
quantum wave packets) with momenta (or momentum- 
space wave packet centroids) p. 

The subtle distinction between 'quantum' momenta q 
and 'classical' momenta p is hardly noticeable and easily 
glossed over in the treatment of a spatially homogeneous 
system like the early universe, but it becomes more ob- 
vious upon consideration of spatial dependence. This is 
because writing down an expression like p(t, x, q) imme- 
diately brings to mind the quantum mechanical incom- 
patibility of position and momentum. In a limit in which 
the neutrinos' motion through spacetime is expected to 
be classical, we want instead an object p(t, x, p) in which 
p represents something other than quantum numbers of a 
momentum eigenstate. If approached from a fully quan- 
tum perspective, obtaining p(t, x, p) will be expected to 
somehow involve a Wigner transformation [l4| (see also 
Ref. [H for a general treatment, having a somewhat dif- 
ferent flavor than that presented in Sec. IIIII below, that 
employs a Wigner transformation in connection with rel- 
ativistic quantum fields). A Wigner transformation is 
a Fourier transformation with respect to a spatial differ- 
ence variable, with p entering as this difference variable's 
'Fourier conjugate.' It is also natural (and correct) to 
guess that a spatial derivative Liouville term p l dpjdx 1 
would appear in the equation of motion for p(t, x, p). 

A few previous efforts towards a kinetic theory of neu- 
trinos with flavor mixing have noted the existence of a 
spatial derivative term in the Liouville operator. These 
might be divided into two broad classes. In several cases 
an explicit derivation of this term is absent [J [H, [H, E3 j 
but its expected presence in a Heisenberg-like equation 
of motion is noted, based on appeals to literature in non- 
relativistic statistical physics [l8l . or quantum optics 
[20| . The other class includes two works [2l|, [HI m which 



3 



spatial derivatives appear in the course of the derivation, 
thanks to use in one way or another of the Dirac equation 
obeyed by neutrino quantum field operators. At least 
from the perspective of working astrophysicists, this sec- 
ond class of approaches has not lent itself to the greatest 
transparency. 

The purpose of this paper is to elucidate the phrase 
"somehow goes over to" in the sentence that follows Eq. 
([S]) , with an emphasis on obtaining the spatial derivative 
and vacuum flavor mixing terms in the flat spacetime 
Liouville equations obeyed by relativistic neutrino and 
antineutrino single-particle distribution matrices (here 
'relativistic' means that only terms of 0(m 2 ,/E v ) are 
kept, where m v and E v are characteristic neutrino mass 
and energy scales). Goals include increased simplicity 
and transparency in comparison with available explicit 
derivations [U and a more detailed discussion of the 
physical interpretation of the Wigner-transformed 'den- 
sity function' (mean value of paired quantum field op- 
erators). The value in this lies in an improved under- 
standing of how classical expressions emerge from quan- 
tum formalisms. This bridge between the classical and 
quantum worlds will facilitate study of the potential im- 
pact of flavor mixing on the decoupling of neutrinos in 
(for instance) core-collapse supernovae: the macroscopic 
quantum effect of flavor mixing must be retained, but 
the neutrinos' motion through spacetime goes over to a 
classical description in order that the system's formal de- 
pendence on time, position, and momentum be simplified 
to a degree comparable to that exhibited by a classical 
single-particle distribution function. (While the formal 
dependence of the neutrino distributions on time, posi- 
tion, and momentum are not unlike those of a classical 
single-particle distribution function, there are indications 
that flavor mixing may lead to new and complicated be- 
havior in the supernova environment [H, [HIH-) 

Two different accounts of the construction of neu- 
trino and antineutrino distribution matrices p(t, x, p) 
and p(t, x, p) and the Liouville equations they obey in 
the absence of interactions are presented in the following 
sections. In Sec. [TT] a simple model of the flavor evo- 
lution of a single neutrino is taken as a starting point. 
In this approach the quantum treatment is restricted to 
flavor evolution; neutrinos are assumed from the outset 
to follow classical trajectories in spacetime. Distribution 
matrices are built up from single-neutrino states, and the 
Liouville equation follows from the evolution equation for 
these individual states. In contrast, a 'density function' 
constructed from quantum fields, and the equations of 
motion it obeys, are the basis of the approach presented 
in Sec. IIIII In this case classical expressions involving 
spacetime variables are not immediately obvious, but can 
be drawn out through use of a Wigner transformation. 
Section IIVI contains a summary and some remarks look- 
ing ahead towards the derivation of interactions from this 

second formalism. Metric signature H and units in 

which h = c = 1 are employed throughout. 



II. STARTING FROM A SIMPLE MODEL OF 
FLAVOR MIXING 

A simple model of neutrino flavor mixing postulates 
the existence of flavor and mass eigenstates. Neutrino 
and antineutrino flavor eigenstates (labeled by a) are re- 
lated to mass eigenstates (labeled by i) by 



i 

\h>w\a) = \vw\i) 



(7) 
(8) 



respectively, where W is the neutrino or antineutrino's 
classical worldline [3£| . These basis transformations fea- 
ture the same unitary matrix that relates neutrino flavor 
and mass quantum field operators: v a (x) = 53* U a i Vi{x). 

Consider a 'Schrodinger picture' in which a neutrino 
state evolves along a worldline W with affine parameter 
A. Let |ia#/(A); a) denote a neutrino that was born in 
flavor a at A = and then translated to A by a uni- 
tary 'worldline evolution operator' ^V(A, 0). As usual 
for unitary transformations parametrized by a continu- 
ous variable, worldline translations are generated by a 
Hermitian operator, here denoted A^y: 



^V(A + d\, A) = 1 -iA^-dA, 
whence the 'Schrodinger equation' 

\vw{\)\a) = Aw \vw(\)\a) . 

A definition of A^- is needed. In flat spacetime 
d . „ d 



(9) 



(10) 



(11) 



where the classical four- momentum py/ is tangent to W . 
The familiar significance of id/dx 11 as a representation 
of the generator of spacetime translations motivates the 
construction of a flavor evolution operator PlL modeled 
on the four-momentum of a particle approaching the rel- 
ativistic limit: 



pO 



W 



P#"| 



M 2 



M 



\Pw\ 



2py 



(12) 
(13) 



where the mass operator M with eigenvalues is not 
diagonal in the flavor basis. Then Ay = p^LJV M , and 
Eq. (HUJ) becomes 

[ 1\ l^( A ) ; a ^ = ~2~ W^W' a ) i (I 4 ) 

a covariant version of the familiar [2~ij | neutrino flavor 
evolution equation. Antineutrino states obey the same 
equation. 
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The off-diagonal terms in the flavor-basis representa- 
tion of M 2 imply that a Schrodinger-picture neutrino 
state | v~g/ (A) ; a) that begins life in flavor a evolves into a 
superposition of all flavors 0. The 'oscillation probabil- 
ity' | {vy) Ww{X)', a) | 2 for a neutrino flavor transforma- 
tion v a — ► up that follows from solution of Eq. (fT4|) is a 
function of A; but it agrees with the usual expression [24[ 
for the vacuum flavor oscillation probability as a function 
of spatial distance L in a frame at rest with respect to 
the source and detector, as can be seen by noting that 
p l w = dx l /dX implies (in flat spacetime) that the world- 
line's affine parameter is equal to Xl = L/\py^\ when the 
neutrino has traveled a spatial distance L. 

The operators giving rise to neutrino and antincu- 
trino distribution matrices can be constructed from states 
\v-g/{X);a) and \v^-(X);a) respectively. The density op- 
erator corresponding to the pure state \vy^(X); a) — that 
is, the operator whose matrix elements comprise the den- 
sity matrix describing a single neutrino with worldline W 
that began life with definite flavor a — is 



so that 



Pw, a (X) = \vw{X);a) {v w (X);a\ 



(15) 



Suppose we have an ensemble of systems of noninteract- 
ing neutrinos with definite flavors a at A = on worldline 
W . Let fw,a(Q) be the ensemble-averaged number of a 
neutrinos at A = on W; then the single-particle neu- 
trino distribution operator describing the ensemble is 



pw (A). 



(16) 



M 2 ap = (v w ;a\M 2 \v w ;0) . 



(20) 



However, to obtain the same representation M^g in the 
antincutrino case the 'backwards' definition 



so that 



M 2 = J2 M *P \ v w,P) {ihr\a\, 

ct,P 



Ml B = (^W,P\M 2 \D^;a) , 



(21) 



(22) 



is required to compensate for the opposite transforma- 
tions of neutrino and antineutrino states in Eqs. ([7]) and 
([8]). If the neutrino and antineutrino distribution opera- 
tors are expanded analogously, 

P*'{X) = ^ Pw,aft{X) \vyr\ a) (i>yr;0\ , (23) 
Pw( x ) = ^2pw, a p(X) |zV;/3) {&w\ot\ , (24) 

a,j3 

then the matrix representations M 2 , p-^-(A), and pyp-(X) 
all transform the same way in species (flavor/mass) 
space: if A represents any of these matrices, then the 
flavor representations are related to the mass represen- 
tations by Afl a vor = U ^4mass U' , where U has elements 
U a i. However, the matrix representations of Eqs. (JTTJ) 
and (|18p now acquire a sign difference: 



Its equation of motion is 

• d - m 1 



M 2 ,p w {X) 



(17) 



which follows directly from Eqs. ([141 and (fT5|) . With 
replacements |z/y(A);a) — > |P*-(A); (A) - 

Pw.aW, /V,«*(0) -> jV, Q (0), and pw(X) -> p w {X) the 
same construction holds, so that 



• d '- m 1 



M 2 ,p^{X) 



(18) 



in the case of antineutrinos as well. 

While the neutrino and antineutrino distribution oper- 
ators obey the same equation of motion, it is convenient 
to define the matrix representations of these equations 
in such a way that that they acquire a relative sign dif- 
ference. The reason is that it is desirable for the matrix 
representations M 2 g of the flavor-basis squared mass op- 
erator M 2 to be the same in the neutrino and antineu- 
trino cases (and equal to the square of the mass matrix 
in the Lagrangian for free neutrino flavor fields). In the 
neutrino case this requirement is consistent with the stan- 
dard construction of a matrix representation: 



i^A~P^( A ) 
i^P^(A) 



1 



M 2 ,p w (X)] , 
-[M 2 ,^(A)] 



(25) 
(26) 



for neutrinos and antineutrinos respectively. Note that 
the absence of hats indicates that these are matrix equa- 
tions rather than operator equations. 

The elements of p-p-(X) and p^r(A) deserve further in- 
spection. They are 

(A) = £)/*7r(0) (v^;a |za^(A);7) 

7 

x(u w (X);j\u^;0), (27) 
Pw,a0{X) = ^ /ir, 7 (0) (i^r; |%'(A); 7) 



x (ia^(A);7 \vw,a) 



(28) 



As the sum of the initial numbers of neutrinos in flavors 
7 on worldline W , weighted by the probabilities of flavor 
transitions 7 — > a, the diagonal elements 

PW ,aa 

(A) = y^\{v w ;a |^(A);7)| 2 /V, 7 (0) (29) 



M 2 = V \ vw - a ) 



a,B 



(19) are equal to fw.a(X), the number of neutrinos of fla- 
vor a at A (and similarly for antineutrinos). Manifestly, 
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the traces of pyp(\) and p^(A) are respectively equal to 
the numbers of neutrinos and antineutrinos of all species 
on worldline W . The off-diagonal elements quantify the 
overlap in flavors a and (3 generated from the initial num- 
bers of neutrinos in flavors 7. 

The desired Liouville equations are close at hand. A 
particular value of A on the worldline W specifies a point 
in spacetime, and the on-shell tangent vector to W coin- 
cides with the neutrino momentum. Therefore, if atten- 
tion is broadened from a single worldline to a collection 
of them forming a congruence of curves in phase space, 
then the specification of W and dependence on A em- 
ployed thus far are equivalent to dependence on i, x, p in 
some coordinate system. Synchronize parametrizations 
in the congruence of curves such that for each worldline 
A = corresponds to t = in a chosen coordinate system; 
then the average particle numbers per worldline /#-, a (0) 
and /ir,a(A) encountered above correspond to / a (0, x, p) 
and f a (t, x, p), where these latter quantities are classical 
distribution functions as in Eq. ([T]). Therefore, with 
a choice of coordinate system, the pw(\) pertaining to 
a set of neighboring worldlines crossing an infinitesimal 
spacelike hypersurface in phase space may be denoted 
p(t, x, p) d 3 x<i 3 p/(27r) 3 , and similarly for antineutrinos. 
(Relativistic neutrinos and antineutrinos produced by 
V — A interactions have spin degeneracy g = 1). 

Hence we have neutrino and antineutrino distribution 
matrices p(t, x, p) and p{t, x, p); taking into account Eqs. 
(fTTjl . (125|) . and |26|) . together with the Liouville theorem 
(invariance of phase space volume elements [6j, LD, |8(), 
these distribution matrices satisfy the Liouville equations 



1 



T 0- 



V 



d 

dx^ 
d 

dx* 1 



p(t,x,p) + i [M 2 ,p(;,x, P )] = 0, (30) 



p(t,x,p)-- [M 2 ,p(i,x,p)] = 0. (31) 



The flavor/mass structure of the Hermitian matrices M 2 , 
p(t,x,p), and p(t, x, p) is given in Eqs. (f2T)|) , and 
([27]) . (|28|) . and their transformation properties are de- 
scribed in the text between these. The diagonal elements 
of the distribution matrices are real and are classical dis- 
tribution functions, as in Eq. fl}, f° r the particle types 
of the chosen representation (flavor or mass). The off- 
diagonal elements quantify the extent of mixing (species 
superpositions) present in the neutrino gas. 

Observable neutrino flavor mixing phenomena depend 
on the differences of squared mass eigenvalues Sji = 
m 2 ~ m 2 , and (aside from upper limits) such differences 
are the only data on neutrino mass that have been ex- 
perimentally determined [2~jj . That flavor mixing prob- 
abilities do not depend on absolute masses (other than 
satisfaction of the relativistic limit) is apparent when the 
squared mass matrix is decomposed as M 2 = £ + A, 
where S is proportional to the identity matrix and A is 
the traceless part. For instance, in the standard case of 
three neutrino species, 



1 



- (m 2 + m\ + ml) 1 (33) 

6 Vo 1/ 



in any basis, and 

(A) mass = (Af 2 ) mass — (£) 



(34) 











I ( -$21 - h\ 
= -[ S 21 -S 32 1(35) 

6 V 632 + S31 ) 

in the mass basis. Because £ cancels out of the commu- 
tator, Eqs. dm)-© become [H 



,1 



dx^ 



p(i,x,p) + -[A,p(i,x,p)] = 0, (36) 
?(i,x,p)-i[A,p(i,x,p)] - 0. (37) 



And going back to Eq. (|T4]) . S merely gives rise to an 
overall phase that cancels in flavor transition probabili- 
ties. 



III. STARTING FROM A QUANTUM 'DENSITY 
FUNCTION' 



A more fundamental approach than that presented in 
the previous section begins with a quantum 'density func- 
tion,' the mean value of a pair of normal-ordered neutrino 
quantum field operators. (Unlike the previous section, 
the convention of denoting operators by hats is aban- 
doned here because distinctions between operators and 
matrix representations thereof using the same symbol 
will not be required.) The focus here is on free neutrino 
fields. Some details regarding these — including conven- 
tions, and reminders about behavior in the relativistic 
limit — are given in the Appendix. 



A. Quantum density function 

Define the free-field quantum 'density function' 
r(y, z) — a function of spacetime positions y and z — in 
terms of the neutrino quantum field operators v(y) and 
D{z): 



(38) 



Tr(M 



f2\ 



(32) 



The subscripts i, j index fields of definite mass, which are 
the 'physical fields' for which the usual quantization in 
terms of Fock states makes sense [32j . The superscripts 
£, m are spinor indices. In this case the bar on v(z) 
denotes the Pauli conjugate v(z) — v*{z)^°; note that 
in most other instances later in this section a bar sim- 
ply labels a quantity related to antineutrinos rather than 
a Pauli conjugate. The angle brackets signify both the 
taking of an expectation value with respect to a many- 
particle quantum state and an average over a statistical 
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ensemble of such quantum states. The ./V on the right- 
hand side denotes 'normal ordering,' which specifies that 
creation operators are to be placed to the left of anni- 
hilation operators, with the introduction of minus signs 
appropriate to the interchange of fermionic operators as 
needed. In particular, separating the neutrino field oper- 
ator 



v{y) = A(y) + B(y) 



(39) 



into its positive- and negative-frequency parts A(y) and 
B(y), the density operator becomes 

irfr(y,z) = - (Af(z)Al(y)) + (Bj(y)B™(z)) . (40) 

As will become clear below, the first and second terms are 
associated with the densities of neutrinos and antineutri- 
nos respectively. Rapidly oscillating cross terms between 
positive- and negative- frequency parts are not relevant to 
the macroscopic limit, and have been dropped (l2l. Il5j. 

In most cases of practical interest the complications 
of spin can be eliminated. This is because the combi- 
nation of V — A neutrino interactions with a relativistic 
limit to 0(m 2 /E v ), where m„ and E v are characteris- 
tic neutrino mass and energy scales, ensures that only 
negative-helicity neutrinos and positive-helicity antineu- 
trinos need be considered. The simplifications that result 
from taking spin transitions off the table are twofold. 

The first simplification resulting from the irrelevance 
of spin pertains to the equations of motion employed. 
While rf™(y, z) obeys the Dirac equation by virtue of its 
construction from Dirac fields, each spinor-space compo- 
nent of this density function also obeys the Klein-Gordon 
equation. Hence, when considerations of spin are irrel- 
evant, the Klein-Gordon equation can be used from the 
outset. At first glance this may seem counterproductive, 
for the same reason the Dirac equation was invented in 
the first place: like Dirac — and with a somewhat related 
motivation, namely the maintenance of positive probabil- 
ity distributions — we are ultimately after equations that 
are first order rather than second order in time. How- 
ever, we shall see that in the present context the desired 
first-order equations emerge very naturally from a combi- 
nation of Klein- Gordon equations, by virtue of a Wigner 
transformation. 

The second simplification resulting from the neglect of 
spin is that the 4x4 spinor structure of Tf" l (y, z) can be 
eliminated, so that focus shifts to entities without spinor 
indices. 



B. Obtaining a first-order equation 

An explicit account of these simplifications begins with 
two Klein-Gordon equations satisfied by the density func- 
tion before the relativistic limit is taken — one with re- 
spect to y, and one with respect to z. Written in matrix 
form with species indices suppressed, these two equations 



□„ T em (y,z) + M 2 T em (y,z) 
n z T lm {y,z)+T lm {y,z)M 2 

where (for instance) 

_d d_ d 2 a 

Dy ~dy~dy^~Wr y 



0. 
0. 



(41) 
(42) 



(43) 



is the d'Alembertian with respect to spacetime position 
y. The difference of Eqs. (J4TJ) and (|42j) is 



(44) 



(□„ ~ □,) r £m (y, z) + [A, T lm (y, z)} = 0, 



in which M 2 has been replaced by its traceless part A 
containing only squared mass differences as described in 
the last paragraph of Sec. |TTJ 

Turn now to the change of variables associated with a 
Wigner transformation. Rewrite y and z in terms of new 
spacetime variables x and 5: 



2' 



x . 

2 



(45) 
(46) 



The meanings of x and S begin to become more apparent 
from the inverse transformations 



2 (» + *)> 



= = y 



(47) 
(48) 



The average spacetime position x will become the 'macro- 
scopic' position variable in classical expressions obtained 
from the present quantum formalism. The difference co- 
ordinate 5 is indirectly related to the 'macroscopic' mo- 
mentum variable in such classical expressions. In partic- 
ular, the Wigner transformation is a Fourier transforma- 
tion with respect to S. The Wigner transformation of 
1^ (y, z) yields the 'mixed representation' Qfj n (x,P) of 
the density function: 



iP-H -p£m 

ij 



(49) 



At this point P is not an on-shell four-momentum. How- 
ever we shall see below that it does basically become the 
'macroscopic' momentum variable in derived classical ex- 
pressions (up to a sign difference between the neutrino 
and antineutrino parts). 

The commutator in Eq. (|4"4"| already looks familiar 
from the Liouville equations obtained at the end of Sec. 
HIl it turns out that the change of variables of Eqs. 
([45]) and (|46|) associated with the Wigner transformation 
starts to bring the differential operator term into more 
familiar form as well. Under this change of variables the 
second-order operator becomes 



□t/ -n. 



8_ d_ 

dE dx' 



(50) 
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which is first order with respect to x. Note also that 



T^(y,z) = 1^(3+ =,*-=) CI, 



Al_ e -^gip(x,P), (52) 



where the transformation in the second line is the inverse 
of that in Eq. (|4"9"|). The mixed representation of the 
density function satisfies 

-2iP»-^Q tm {x,P)+ [A,g* m (x,P)] = 0, (53) 

which follows from substitution of Eqs. ([50)) and (|52")) 
into Eq. (f44|) . 



C. Interpretation of the mixed representation 

Consider next the spinor-space structure of the density 
function in the case of practical interest described above, 
in which the neutrino and antineutrino populations are 
overwhelmingly relativistic. In this case it follows from 
the explicit expressions in the Appendix that the spinor 
space structure of the density function reduces to 



^ m {y,z)) 



T LR (y,z) 




(54) 



Note that in the relativistic limit only one 2x2 block is 
nonzero. The notation T LR (y, z) for this block denotes 
the fact that it would be the block projected out if T(y, z) 
were sandwiched between the left- and right-projection 
matrices P L and P R of Eqs. (|A.7|) and (|A.8|) . Explicitly 



^(y,z) 



f d 3 g d 3 u r i(uj . z _ qi . y) lr, , 
J (2tt) 3 (2tt) 3 L iN « lq ' Uj 

(55) 



(2tt) 3 (2tt)' 
e H«-J/-«i^)N§ R (q, u ) 



Here (<#) = (£<,,*, q) with P q ,, ; = y/]q\ 2 +mf « |q| + 
m?/2|q| (and similarly for the components of Uj), and to 
0(ml/E v ) 



N 



Li? 



(q,u) = {ai tU a n , U i), 



(56) 

Nf/(q,u) = vWJ (bl u bu,r,j), (57) 

as discussed in the Appendix. (Note that £^ and r/^ are 
two-component spinors.) 

Moving to the mixed representation provides a con- 
venient means of separating the neutrino and antineu- 
trino parts of the density function, by making manifest 
its positive- and negative- frequency parts. Apply the 
Wigner transformation of Eq. (|49| to Eq. (|55|) and find 



d 3 q d 3 u 

(2tt) 3 (2tt) 3 



e -i(«- UJ ).x N £fl( qi u ) 



x (2^) 4 <5 4 (P - 2LJ 



x (2^) 4 <P P + 



<? 4 + % 



(58) 



Because g.; and itj are on-shell momenta, it is evident that 
P° > in the first term and P° < in the second term. 
Separate neutrino and antineutrino density functions are 
obtained by projecting out these positive- and negative- 
frequency parts with step functions 8 (P°) and (— P°) : 



iG% R {x,p) 



iJd 4 P6\ P -P)8(P°) gf R {x,P 
d 3 q d 3 u 



(2tt) 3 (2tt) 
x (27r) 4 ^ 4 



3 e- i ( 9 *-^) ,a; N^ ii (q,u) 



(59) 



and 



iGf R (x,p) 



i J d 4 P <5 4 ( P + P)6 (-P ) gjf-{x, P) 



x (2tt) 4 S 4 ( p - 



(60) 



Note that delta functions included in the projection op- 
erations yield the momentum label changes P — » p in 
the case of the neutrino density function G^ R {x,p) and 
P — > — p in the case of the antineutrino density function 
Gf R (x,p). These density functions satisfy 

-2ip^G LR (x,p) + [A,G LR (x,p)] = 0, (61) 

2ip^G LR (x,p) + [A,G LR (x,p)] = 0, (62) 

which follow from applying to Eq. (|53p the same projec- 
tions appearing in Eqs. (|59|) and ([60]) . 

In general cases the mixed representation provides 
complementary position and momentum probability dis- 
tributions, as required by the quantum mechanical in- 
compatibility of position and momentum, rather than a 
joint position/momentum distribution. 



1. Position distribution 

A position distribution is obtained by integrating the 
mixed representation over all 'momenta.' This can be 
seen by comparison of the diagonal elements of Eqs. (|59"|) 
and (f6"0"|) in species space with the relevant component of 
the (normal-ordered) number current (Nvi(x)^Vi{x)) of 
species i. 

In particular, the net neutrino number density of 
species i — that is, the difference between the neutrino 
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and antineutrino number densities or position distribu- 
tions rii(x) and fii(x) — is the Oth spacetime component 
of the number current: 

rii(x) - ni(x) = {Nvi{x)^°Vi{x)) . (63) 

More explicitly, 

d 3 q d 3 u - 



rii(x) - rii(x) 



(2tt) 3 (2tt) 3 



-ife-«*)-^(q !U ) 



'^^(q.u) 



where 



M(q, u) 
M(<1, u) 



(64) 

(65) 
(66) 



Aside from the single species index i compared with the 
potentially different indices i and j, the difference be- 
tween Eqs. (|65|) . ((66)) and (|56|l . ((57)1 is that the one pair 
has inner products of two-component spinors, while the 
other has an outer product giving rise to a 2 x 2 matrix 
in spinor space. 

The neutrino and antineutrino contributions to Eq. 
are readily obtained from the species-space diago- 
nal components of Eqs. ((59)) and |60)) respectively. Inte- 
grating over p, and using the fact that the inner product 
of any two spinors is equal to the trace of their outer 
product, one finds that 



n,i{x)-m(x) = - 



d 4 p 



TT[iG^(x,p)+iG^ R (x,p)] , 

(67) 



where the trace is over the spinor indices of the 2x2 
blocks. This motivates the definition of spatial 'number 
density matrices' 



/(04 1V [- i ^(*,--P)] 



(68) 



and 



^ Tr[iGf/(i,x,p)], (69) 



whose diagonal elements are spatial number densities 
of the various neutrino and antineutrino species re- 
spectively. (The dependence on spacetime components 
(x M ) = (t, x) has been displayed here more explicitly.) 



2. Momentum distribution 

A momentum distribution is obtained by integrating 
the mixed representation over the volume of the system. 
Integrating Eq. (f6~4")) over all space results in a factor 



(27r) 3 (5 3 (u — q), so that the net total neutrino number of 
species i is 

Ni-Ni = / d 3 x [m(x) - fiiix)] 



(27T) 



(70) 



Comparison with the definition of an occupation num- 
ber in Eq. ([5]) — whose integral over q also gives a total 
number of particles in the system — implies that the dif- 
ference of neutrino and antineutrino occupation numbers 
or momentum distributions is 

m{t,q) -fii(t,q) = ^ ((a q .u a q,l,*) ~ (&q,t ,i b i,T ,*)) ■ 

(71) 

Precisely the same expression, but with momentum label 
q replaced by p, is obtained from the diagonal elements of 
Eqs. ((59)) and ((60)) by integrating over x, integrating over 
p° (which puts the momentum p on shell), and making 
use of Eqs. (|A.29p and (|A.30[) in the Appendix: 



n,(t,p) - rii(t,p) 



V 



dp" 
W) 



Tr[iG^ R (x,p)+iG^ R (x,p)} . 

(72) 



This motivates the definition of 'occupation matrices' 



Pij 
and 

Pij 



(t,p) = 1 Jd^J^ Tr [-iGffit^p)} (73) 
(t,p) = IJd^J^ Tr [iG^(t,x,p)] , (74) 



whose diagonal elements are occupation numbers of the 
various neutrino and antineutrino species respectively. 
(Here the components of p = q = u are quantum num- 
bers of momentum eigenstates, a role denoted by q in 
Sec. U) 



3. Classical joint position/momentum distribution 

Beyond the complementary position and momentum 
probability distributions generally available from the 
mixed representation of a density function, a joint posi- 
tion/momentum probability distribution — akin to a clas- 
sical one-particle distribution function or phase space 
density — is expected to emerge when the spacetime and 
momentum dependences of a neutrino ensemble satisfy 
appropriate conditions. In particular a classical system 
is characterized by 



£T»1, PL>1, 



(75) 
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where E and P are characteristic energy and momentum 
scales and T and L are characteristic time and length 
scales. 

In order to elucidate this classical limit it is neces- 
sary to examine N^- fl (q, u) and N^- fl (q, u) in Eqs. ([59]) 
and ([60]) more closely, beginning with a specification of 
the state with respect to which the expectation values in 
Eqs. and ([57)) are taken. For illustrative purposes 

a pure state |$jv) 01 N neutrinos will be discussed here. 
The extension to systems represented by pure or mixed 
states which also include antincutrinos or other particle 
types might provoke complications (or at least changes) 
in notation but would involve no significant additional 
conceptual difficulties. 

A pure state |$jv) of N neutrinos is built up out 
of multiparticle momentum eigenstates |kiii . . .h^i^), 
where the k are momenta and the i label mass eigenval- 
ues. These momentum eigenstates result from the action 
of antisymmetrized products of neutrino creation opera- 
tors a k , i upon the vacuum, together with energy factors 

\/2Ek, such that these antisymmetric states are normal- 
ized according to 

(kiij . . . k' N ,i' N ,\kiii . . . k N i N ) 

N 

= ^^^n^kJ^Tr)^ 3 ^ - k Q )^ a . ia .(76) 

V a=l 

The sum is over all permutations V of the list of parti- 
cle labels 1 . . . N indexed by a, with S-p equal to 1 for 
even permutations and —1 for odd permutations, while 
Va is the particle label moved to the ath position un- 
der a particular permutation V . This sum reflects the 
complete antisymmetry of these multi-neutrino states in 
that it vanishes if any two momentum/mass label pairs 
are equal. 

In accordance with the assumption that the neutrino 
gas can be represented by a density function — that is, 
that correlations can be ignored — take |$jv) to be a su- 
perposition of momentum eigenstates constructed with 
N independent single-particle wave packets: 

x |kiii . . . kjvijv) • (77) 

The </> p (k) are real- valued wave packet 'envelopes' cen- 
tered on p, each normalized such that 

r ^3), 

J [0p(k)]2 = L (78) 

According to the usual wave packet technology these 
momentum-space envelopes are related to position-space 
wave packet envelopes V>x , P (x), peaked about x , by 

/r/ 3 k 
^0 p (k)e ik -( x --), (79) 



whose normalization 

|d 3 x^ X0 , P (x)| 2 = l (80) 

follows from Eq. ([78]) . (For instance, if the </> p (k) are 
taken to be Gaussian, then 

^ X0 , P (x)=^ XD (x)e i P-( x - x °), (81) 

where ^ XD ( x ) is a (real- valued) Gaussian centered on xo.) 

Before explaining the particular choice of superposi- 
tion of mass eigenstates represented by the sum over 
i a in Eq. |77|) it is helpful to clarify that |$jv) is a 
Heisenberg-picture state which for dcfinitcncss is taken 
to correspond to a collection of N relativistic neutrinos 
at t = 0. These neutrinos are somewhat localized in both 
momentum space and position space, the ath neutrino 
being localized around p a and xo^a 

respectively (where 
the subscript is a reminder that this position localiza- 
tion is that which applies at t = 0). Hence the particular 
superposition ^a a i a applies to a system in which the 
ath neutrino is in 'definite flavor' a a at t = 0. This corre- 
sponds to the assumption used for the sake of illustration 
in Sec. |TI]that all the neutrinos were in definite flavors a 
at t — 0. More generally the various neutrinos could be 
taken to be have been created in definite flavors (that is, 
in association with different charged leptons) at various 
different times t a < 0; in this case a more general su- 
perposition J2i c a a i a with J2i \ c a a i a \ 2 = 1 would apply, 
where c aa i a encodes the amplitude for the ath neutrino to 
be in mass eigenstate i a after evolution from its creation 
in flavor a a at t a < until t = 0. 

In order to obtain expectation values with respect to 
|$iv) it is necessary to know its norm, and this in turn 
requires an understanding of how the Pauli exclusion 
principle — manifest in the complete antisymmetry of the 
multiparticle momentum eigenstates |kiii . . . kjviiv) — 
impacts our localized particles. Consider a situation in 
which (at least) two neutrinos have the same flavor con- 
tent, and also wave packet envelopes </> p (k) in Eq. ([77]) 
with identical shapes, centroids p, and spatial offsets 
x . In this case |$jv) vanishes, because a product of 
wave packet envelopes/mass amplitudes that is symmet- 
ric with respect to (at least) two sets of quantum num- 
bers k, i is 'contracted' with the completely antisymmet- 
ric state |kiii . . .kjvijy)' More generally, the norm of 
|$jv) is 

($Ar|$Ar) = 1 

N [ d 3 k a 

a=l J ^ 

Xe ik.'(x P .,o-x., ) ) ( 82 ) 

which follows from Eqs. ([76]) - (|78p and the unitarity of 
the mixing matrix. The sum is now over all permuta- 
tions except the identity. In the aforementioned case of 
perfect wave packet overlap this sum becomes — 1, and 
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($Ar|<I>Ar) vanishes. At the other extreme, if there is ab- 
solutely no wave packet overlap in either position or mo- 
mentum space, then the sum in second and third lines 
of Eq. (|82[) vanishes, so that ($^1$^) is unity. Between 
these extremes partial wave packet overlaps give this sum 
a value somewhere between and —1, and in turn the 
norm takes continuous values between 1 and 0. There- 
fore 'wave packet smearing' has a graduated impact upon 
manifestations of the exclusion principle in 'microscopic' 
representations of localized particles. 

However, with an end goal of a macroscopic treat- 
ment of a neutrino gas it is neither possible nor desir- 
able to follow the details of all particles' wave packet 
shapes and overlaps; what is of interest instead is how the 
Pauli exclusion principle percolates down to the classical 
limit. To this end take a view that is sufficiently 'coarse- 
grained' as to impose a binary distinction between com- 
plete overlap and complete non-overlap of momentum- 
and position-space wave packet envelopes; then Eq. ([52|) 
goes over to 

N 

(2^) 3 <5 3 ( 

P"Pa Pa ) 

7V1 a=l 

x S 3 (x Vafl - x a , ). (83) 

Here the normalizations associated with the momentum 
and position 6 functions mirror Eqs. ([75)) and ([50]) re- 
spectively; in particular the S functions of zero argument 
are to be interpreted as (2n) 3 S 3 (p-p a — Pa)| PPa = Pa = V 
and <5 3 (x-p a ,o - x Qi0 )|x T , Qi0 =x a , = V~ l , where V is the 
(effectively infinite) quantization volume. Hence the ex- 
clusion principle is promoted from strict applicability to 
global momentum eigenstates to 'for all practical pur- 
poses' joint applicability to the centroids of momentum- 
and position-space wave packets [4lJ. (Harking back to 
the paragraph before last, the factor 8 ava a a pertains to 
the special initial condition in which all neutrinos have 
definite flavor at t = 0; more generally, this Kronecker S 
would be replaced by J2i a c a Va i a c a a ia-) 

With this sort of multiparticle state in mind the expec- 
tation values in Eqs. (|56|) and (|57|) can be evaluated. In 
particular, continuing with the example of a pure state 
of N neutrinos given by Eq. (|77|) , the action of an anni- 
hilation operator upon |$jv) results in a sum of N terms: 

N 



where 



(84) 



N 

n 



d 3 k b Mkt) -tk^,, 



x |ki«i . . . ko-ita-i, k a+ ii Q+ i . . . k N i N ) (85) 

is the N — 1 particle state resulting from the 'removal' of 
the ath neutrino from The sum in Eq. (|84p arises 



because of the anticommutation rule obeyed by neutrino 
creation and annihilation operators, together with the 
fact (mentioned above) that the multiparticle momen- 
tum eigenstates |kiii . . .h^i^} are constructed by act- 
ing upon the vacuum with antisymmetrized products of 
neutrino creation operators ^ v (In particular, the re- 
lation 

a q ,l,i | kj.ii • ■ ■ k/vijv) 

N 

= ^(-) a+1 (27r) 3 v^i; ( 5 3 (q-k a ) ( 5 lla 

a=J 

x |kiii . . . k _ii _i, k a+ ii a+ i . . . k N i N ) (86) 



has been employed in obtaining Eqs. ([84]) and l|85p.) 
The expectation value in Eq. (|56p . which follows from 
Eq. ((SI]), is 



N N 



(<u a ^) = EE(-) b+1 (-) a+1 v( u )^(q) 



6=1 a=l 

x e iu Xi, -° i 



-iq-x a , tj . rr* 



x (® Nji \<S> N( i) / (®n\<5>n) 



(87) 



Evaluate the norms ($^1$^) and (^jvl'&Ar) in the 
coarse-grained 'all-or-nothing' view of wave packet over- 
lap that led to Eq. ([53")) , and let all the neutrino fla- 
vors and wave packet momentum centroids and positions 
be non-overlapping as required for (f&jvl'&jv) — > 1 in- 
stead of ($jv|$jv) - * 0. In this same approximation 
(*jvjil*jVf() -> S ab . Hence Eq. flSTJ becomes 



( a u,i,j a q,bi> 



N 

E 

a=l 

x e 



p o (u)0 Pa (q) 

m-x..o e -iq-x..o n {J* 



when the norms are evaluated in the coarse-grained posi- 
tion/momentum picture, with all of the momentum and 
position wave packet centroids p and Xo being different 
by at least a wave packet width or so. 

Before applying the coarse-grained 'all-or-nothing' 
view of wave packet overlap to the remaining wave pack- 
ets in Eq. ([55)) it is appropriate to consider the evolution 
of these wave packets that follows from putting this ex- 
pression for (ajj | jO-<i,l,i) back into the density function's 
mixed representation of Eq. ([5"9"]) . via Eq. ([5S|) . In this 
connection it is convenient to 'open up' the 5 function in 
Eq. (J5T 



(89) 



According to the usual wave packet technology, the inte- 
gral over q in each of the N terms in Eq. ([56]) yields a 
moving position-space wave packet: 



^L^ Pa (q) e-^-^De-^-^ 
' ^x a(t+S o/ 2 )- 3 / 2 (x) e- ip ^< x+ 9) e" 



-ip o -X o , ti 
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Here = (£ Paji ,p a ), and 



Xa(t) = X Q:0 + V a t, 



(91) 



where v a = p a /E Pat i w Po/|Po| is the wave packet's 
group velocity, approximated to 0(rn~/ E v ). The wave 
packet centroid does not follow this classical trajectory 
for generic values of S (all of which are probed, thanks 
to the integral over H in Eq. (J55])); instead, the nota- 
tion ^/ , x a (t+H°/2)+H/2( a; ) f° r the wave packet envelope is 
meant to convey the fact that the centroid follows the 
non-classical trajectory 



7^0 

2 Xa -°-2 +Va { t+ Y 



(92) 



Similarly, the integral over u in each of the N terms in 
Eq. yields a factor 



d 3 u 



(2tt) 3 Vp » 

« C o(i -H0/ 2 ) + B/2(^ e ip --(^t) e'P— 4 f (93) 

which has three differences from Eq. (|90|) : complex con- 
jugation, the replacement S — ► — H, and mass index j 
instead of i. 

Now that the evolution of the wave packets has been 
considered the full expression for the mixed representa- 
tion of the density function can be evaluated. 

The first step is to apply the classicality conditions of 
Eq. (f75|) . which can be related to assumptions of 'slow 
change' and 'weak inhomogeneity. ' Use of Eqs. 
I (ESI), and (J93|) in Eq. (59|) yields 



Pa,i+Pa 



d 4 E D a (x)e^ p '^^)- S 



N 

E 

a=l 



x e- 1 ^- 1 -^*'^* U aaj U* 



(94) 



where 



D a {x) = ^ a(t _ S o /2)+s/2 (x) i> Xa ( t +S°/2)-B/2{%)- (95) 

Because of the wave packets' localization, D a (x) peaks 
at 5 = 0, for which the centroids of both wave packet 
envelopes follow the classical trajectory of Eq. (l9"T]) . Ex- 
panding about 5 = 0, 



D a (x) = |^x (t)(z)| + s 



(96) 



Note that the first correction term, in combination with 
the first exponential in Eq. (|M|) . can be expressed 



Tx D ^ x) 



<9p 



H=0 



i(p- 



.(97) 



In a collection of particles satisfying the classicality con- 
ditions of Eq. (|75|) . it can be expected on dimensional 
grounds that this and higher-order corrections can be ne- 
glected in the sum over a large number of particles N in 
Eq. ([94]) , thanks to the combination of derivatives in Eq. 
(HZ]). 

A few more simple steps bring the position and mo- 
mentum dependence of Eq. (|94p into fully classical form. 
Keeping only the first term in Eq. (|96|) , the integral over 
S in Eq. (|94p yields a four-momentum 8 function (and 
integration over p° reduces this to a three-momentum S 
function). Moreover, in the coarse-grained view of 'all- 
or-nothing' wave packet overlap, this first term of D a (x) 
approaches a sharp restriction to the classical trajectory 
ofEq. JST): 



D a (x) » \ip Xa ( t )(x)\ -» <5 3 (x-x a (t)). 



(98) 



Finally, because the two spinors in the outer product 

£pa£pa P er ta m to the same momentum, Eq. (|A.25[) ap- 
plies; hence taking the trace over this 2x2 block in spinor 
index space results in a factor of unity. All together, the 
neutrino distribution matrix obtained from Eq. (|94p is 



(i,x,p) = ^ x,p)] (99) 

N 

- £(27r) 3 <5 3 (x-x (*))«5 3 (p-p a ) 
l 

e -i(^P..i-.Bp..i)* U aaj U* a t . (100) 



X 



A similar argument applies to antineutrinos, with the 
result 

/3tf(*,x,p) = |^Tr[-iGf/(t,x,p)] (101) 

AT 

-> 5>7r) 3 <5 3 (x-x a (t))<5 3 (p-p a ) 

a=l 

-Sp«,i)* ^ j;* (102) 



for the antineutrino distribution matrix. Here p takes 
on the values of momentum-space wave packet centroids 
(which ultimately correspond to the momenta of classical 
particles), in accordance with the role denoted by p in 
Sec. |H 

The position and momentum dependence of Eqs. (|100|) 
and (|102[) is precisely that of a collection of free classi- 
cal particles, and the associated species-space structure 
is just as expected. These expressions were derived in 
terms of mass fields, because these represent the 'physi- 
cal particles' for which creation and annihilation opera- 
tors obeying appropriate anticommutation relations exist 
for arbitrary momentum. But once an average is taken 
for a system containing only relativistic neutrinos, the 
notion of a 'flavor basis' for the distribution matrices 
becomes workable [42]. From Eq. it is apparent 

that density function Y£k{y,z) constructed from flavor 
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fields v a (y) = U a t vi (y) is related to that constructed 
from mass fields by 



r^(y,z) = J2u ai rf™( y ,z)u 



fir 



(103) 



1,3 



This propagates down to a 'flavor basis' representation 
of the distribution matrices: 



p a p(t,x,p) = ^2U ai p ij (t,x,p)U^ j , (104) 

i ,3 

p af3 (t,x,p) = ^2u ai p i:j (t,yi,p)U^ j . (105) 

i,3 

In particular the diagonal elements can be expressed 

N 

P W («,x,p) - ^^(2^) 3 <5 3 (x-x a (i))<5 3 (p-p Q ) 



a— 1 i 



X 
JV 



Uf3i exp -i 



"2|p.| 



/ 



, (106) 



fcj/jftx.p) -> ^^(2^) 3 (S 3 (x-x a (i))^(p-p a ) 



C^li exp ( -i 



"2Tp7i 



(107) 



These expressions describe collections of neutrinos and 
antineutrinos that begin in flavors a a at t = and fol- 
low classical spacetime trajectories, with the expectation 
values of their flavors along those trajectories varyin g ac - 
cording to familiar vacuum oscillation probabilities [24| . 

To conclude this subsubsection, it is appropriate to 
remark on the illustrative character of this derivation 
of the classical limit of the position/momentum depen- 
dence. A single pure state of neutrinos has been singled 
out for detailed discussion, and the average denoted by 
angled brackets in Eq. (|56|) has been interpreted as an ex- 
pectation value with respect to this pure neutrino state. 
Hence the derivation of Eqs. (|100[) and (|102[) constitutes 
a demonstration of how a single 'microstate' of uncorre- 
cted single-particle quantum wave packets corresponds 
to a single 'microstate' of definite classical trajectories 
(while retaining the quantum mechanical phenomena of 
flavor mixing and Fermi statistics). This adequately il- 
lustrates the issues involved in a passage to the classi- 
cal limit. But it should also be noted that in order for 
p(t, x, p) and p(t, x, p) to actually be single-particle dis- 
tribution matrices for the 'macrostate' of the gas, the 
average represented by angle brackets must be promoted 
to an average over a mixed state of neutrinos or antineu- 
trinos respectively, together with an ensemble average 
over an appropriate statistical distribution of these mixed 
states. (The mixed states of neutrinos or antineutri- 
nos, represented by density matrices spanning all possible 
numbers of neutrinos or antineutrinos with all possible 
independent single-particle wave packets, would be ob- 
tained by integration of the pure states of the entire sys- 
tem over the degrees of freedom of all particle types other 



than neutrinos or antineutrinos, in accordance with the 
standard meaning of mixed states and density matrices 
[25|,1 Then, in an appropriate limit along the lines pre- 
sented here, the position/momentum dependence of the 
single-particle distribution matrices would correspond 
to the position/momentum dependence of single-particle 
distribution functions for a statistical ensemble of classi- 
cal particles. Finally, as a last reminder of another illus- 
trative feature, the factor U aa j U* ai in Eqs. (jlOOp and 
(|102[) corresponds to the special case in which the neu- 
trinos or antineutrinos respectively have definite flavors 
at t = 0. In more general cases one would have c*^ c aa i, 
where the meaning of the amplitude c aa i was explained 
in the fifth paragraph of this subsubsection. 



D. Liouville equations 

The Liouvillc equations for the neutrino and antineu- 
trino distribution matrices of Eqs. (jlOip ('mass ba- 

sis') and (|104[) . (| L 05[) ('flavor basis') follow immediately 
from Eqs. (|6ip and (|62p . In matrix form, with species 
indices suppressed, the Liouville equations are those ob- 
tained at the end of Sec. [Til 



dx^ 
dx^ 



p(t,x, p) + - [A,p(i,x, p) 
p(t,x, p) - - [A,p(t,x, p) 



0, (108) 
0. (109) 



As mentioned immediately following Eq. (|91j) . the tra- 
jectories x a (t) appearing in Eqs. (|100[) and p02[) are null 
to 0[m} v j Ey). In accordance with this — and in order for 
these explicit expressions for p(t 7 x, p) and p(t, x, p) to 
satisfy Eqs. (| 1 08[) and (|109[) — the momentum p in the 
first term of the Liouville equations above should be ap- 
proximated as = (|p|,p). Only in phases in the ex- 
plicit expressions for p(t, x, p) and p(t, x, p) encountered 
in the previous subsection are corrections of C(m 2 /2|p|) 
to E p retained. 



IV. CONCLUSION 

The calculation of neutrino decoupling from dense nu- 
clear matter requires a transport formalism capable of 
handling both collisions and flavor mixing; and the first 
steps towards such a formalism are the construction of 
neutrino and antineutrino 'distribution matrices,' and a 
determination of the Liouville equations they satisfy in 
the noninteracting case. These initial steps have been 
accomplished in two new ways in this paper. Both ap- 
proaches arrive at neutrino and antineutrino distribution 
matrices p{t, x, p) and p(t, x, p) whose dependence on 
time t, position x, and momentum p is classical. In- 
deed the diagonal elements of these distribution matrices 
are classical distribution functions, in the sense of Eq. 
(HJ), for the various neutrino species. The off-diagonal 
elements encode information on species overlap in the 
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neutrino ensemble. The flat-spacetime Liouville equa- 
tions satisfied by p(t, x, p) and p(t, x, p) are given both in 
Eqs. p6 ]) . ([37 ]l and (fTU8|) . (fT09l) . In addition to the usual 
flat-spacetime directional derivative along the phase flow 
p^d/dx* 1 (see also endnote [4Cj), the Liouville operators 
for neutrinos and antineutrinos with flavor mixing in- 
clude a commutator with a matrix containing differences 
of squared neutrino masses (see Eq. (|35p ). with a differ- 
ence in sign between the neutrino and antineutrino cases. 

In the approach of Sec. [IT] the neutrino positions and 
momenta are taken to be classical from the outset: only 
neutrino mass/flavor are treated quantum mechanically. 
Distribution matrices are constructed from the states of 
a covariant version of the familiar simple model of flavor 
mixing, and the Liouville equations follow straightfor- 
wardly from the 'Schrodinger equations' describing the 
evolution of flavor along classical worldlines. 

The second approach — presented in Sec. [TTT] — employs 
a 'density function,' the mean value of paired neutrino 
quantum field operators (Eq. ([38)) ); therefore the classi- 
cal position/momentum dependence must be derived as 
a limit. The key to this is the 'mixed representation' of 
the density function obtained by a Wigner transforma- 
tion (Eq. (|49|) ). By definition the spacetime variable x of 
the mixed representation is the average of the field oper- 
ators' position variables (Eq. (jTf)) ). and the momentum 
variable P of the mixed representation also turns out to 
be (up to a sign in the case of antineutrinos) an average 
of the momenta appearing the field operators' plane wave 
expansions (note the 8 functions in Eqs. and |M|) V 
In order to make a concrete physical connection between 
these density functions and distribution matrices, con- 
siderable space is devoted in Subsec. IIII CI to explicit 
examination of the general relationship of the mixed rep- 
resentation to complementary space and momentum dis- 
tributions, and especially to a detailed illustration of how 
a suitable state of uncorrelated wave packets, satisfying 
basic classicality conditions of slow variation and weak 
inhomogeneity (Eqs. (|75p ). corresponds to a microstate 
of a gas whose neutrinos follow classical trajectories while 
exhibiting the usual flavor oscillations. This last illustra- 
tion is by far the paper's densest thicket; by comparison, 
almost magically simple and elegant is the emergence in 
Subsec. IIIIBI of the Liouville operator from a Wigner 
transformation of the difference of Klein-Gordon equa- 
tions obeyed by the density functions. 

Given the existence of the simpler approach of Sec. 
ITU the question arises as to why it is worth bothering 
with the more fundamental approach of Sec. IIII1 The 
reason is that it is necessary to go beyond the case of 
noninteracting neutrino distributions that satisfy Liou- 
ville equations: neutrino interactions need to be consid- 
ered, and our understanding of neutrino interactions is 
based on quantum field theory. It is true that one can 
attempt to compute neutrino effective masses and inter- 
action rates independently and then insert them into a 
formalism based on the simple model of neutrino flavor 
mixing. But there is a history of overlooking important 



aspects of the problem when such an approach is taken. 
For instance, early works that considered contributions 
of neutrino-neutrino forward scattering to effective neu- 
trino (squared) masses, such as Refs. J2(| H3|, failed 
to recognize the existence of off-diagonal contributions 
[HI, [5^1 . Another example is an apparent failure [l7| to 
recognize that the placement of neutrino blocking factors 
in neutrino interaction rates is nontrivial, due to the fact 
that neutrino distributions are now represented by non- 
commuting matrices (l2j . Such issues are automatically 
raised and naturally handled in a formalism that begins 
by treating all aspects of the problem in terms of quan- 
tum field theory from the beginning; hence, in addition 
to being conceptually satisfying, an approach of this kind 
is also theoretically safe. 

It may then further be asked why a new treatment 
might be desirable if in fact interactions have already 
been responsibly addressed in the literature [l2j]. While 
the handling of general classes of interactions is outlined 
in Ref. fl2l |. the particular interactions relevant to neu- 
trino decoupling are not spelled out with the degree of 
explicit specificity needed by those developing large-scale 
simulations involving neutrino transport. In addition to 
the present work's goal of more thorough understandings 
of the spatial derivative term in the Liouville equation 
and the nature of neutrino distribution matrices, it also 
serves as a first step towards the handling of neutrino 
interactions with a diagrammatic approach based on a 
nonequilibrium Green's function [3(| (see also Ref. [22|). 
In addition to a Green's function, this approach — which 
is sometimes called 'Keldysh theory' — involves the den- 
sity function considered in this paper and two other types 
of field operator pairings. Because it is a diagrammatic 
approach while that of Ref. [l2j is purely algebraic, it 
can be expected to simplify the task of explicitly work- 
ing out the full panoply of interactions needed by those 
wishing to incorporate neutrino flavor mixing into neu- 
trino transport computations. This will be pursued in a 
separate work. 



APPENDIX 

The approach in Sec. IHII to neutrino distribution ma- 
trices and the Liouville equations they obey is rooted in 
quantum field theory. Quantum 'density functions' con- 
structed from neutrino fields simplify considerably in the 
relativistic limit of practical interest. 

Neutrino interactions respect lepton number in the 
Standard Model, but with good evidence from flavor mix- 
ing observations and experiments that neutrinos have 
mass, it is not clear that neutrinos actually carry lep- 
ton number (24|. Standard Model neutrino interactions 
are of the V — A form and therefore involve only left- 
handed neutrino fields. For massless neutrinos this im- 
plies the existence of only two neutrino states for each 
flavor: negative-helicity neutrinos and positive-helicity 
antineutrinos. In the relativistic limit of extensions that 
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give the neutrinos mass, these two states could be either 
left-handed 'neutrino' and right-handed 'antineutrino' 
states associated with a Dirac field, or the left- and right- 
handed states of a self-conjugate Majorana field carrying 
no net lepton number. Differences between these two pos- 
sibilities occur only at 0(rr? v j E 2 .), where m„ and E v are 
characteristic neutrino mass and energy scales. Systems 
in which effects at this scale are relevant are not con- 
sidered in this paper; instead, only terms of 0{m 2 /E v ) 
are kept, which capture the flavor mixing physics rele- 
vant to the neutrino decoupling problem. For definitc- 
ness the language and formalism of Dirac fields are used 
here — that is, 'neutrinos' and 'antineutrinos' and associ- 
ated operators a qir , b^ r and so forth are spoken of — but 
the results to 0(m%/E v ) are the same as if Majorana 
fields were employed. 

The operator ip (y) representing a noninteracting spin- 
1/2 field possessing one or more quantum numbers dis- 
tinguishing particles from antiparticles — that is, a Dirac 
field — is 



^(y) = A i (y) + B e (y), 



(A.l) 



where A e (y) and B e (y) are functions of spacctime posi- 
tion y and constitute 'positive-' and 'negative-frequency' 
parts respectively: 




XX fa'' 1 

r 



^a q , r ,(A.2) 
(A.3) 



Here £ is a spinor index and r labels spin states. The 
momentum 4- vector q has components (q' 1 ) = (i? q ,q), 
where E^ = y/\q\ 2 + m 2 is the on-shell energy of a 
particle of mass m. The positive-frequency term con- 
tains particle annihilation operators a^ r and momentum- 
space Dirac spinors z/(q, r), and the negative- frequency 
term contains antiparticle creation operators 6 q r and 
momentum-space Dirac spinors i/(q, r). The free field 
satisfies the Dirac equation 



if 



d 
dyf" 



A ip{y)- 



(A.4) 



The Dirac spinor indices on 7 M and ip{y) have been sup- 
pressed. Here the Dirac matrices 7 P satisfy the an- 
ticommutation relations {7^,7"} = 2r) fiv , where if ,v 
is the Lorentz metric. The conventions of Ref. [31| 
are followed for units (H = c — 1); metric signature 

(H ); creation/annihilation operator anticommuta- 

tion relations [{a q ,r,4 jS } = {frq,r,&u. s } = (27r) 3 <5 3 (q — 
u)<5 rs , with all other anticommutators vanishing]; single- 
particle states (|q, r) = y / 2~E q "a q r |0)) and their normal- 
ization [(q,r|u,s) = 2_E q (27r) 3 5 3 (q — u)# rs ]; and Dirac 
matrices, which have the 2x2 block form 



(V) 







a' 1 




(A.5) 



Here (cr M ) = (l,c) and (cr M ) = (1, — er), where cr are the 
standard 2x2 Pauli matrices. Moreover the matrix 



(7 5 



-1 



appears in the left and right projection operators 
Pl = ^(l-7 5 ), 



Pr 



;(l + 7 5 ), 



(A.6) 

(A.7) 
(A.8) 



which select the upper two and lower two components of 
Dirac spinors respectively. 

The simplifications incident to the relativistic limit are 
manifest in explicit expressions for the neutrino and an- 
tineutrino momentum-space spinors. These can be writ- 
ten in block form as 



v(q,r) 



(A.9) 
(A.10) 



It is convenient to define the two-component spinors £ q 
and ?7 q in terms of eigenspinors x q that satisfy the rela- 
tions 



(A.ll) 



The Dirac spinors are associated with left-handed 
(negative-helicity J,) and right-handed (positive- hclicity 
t) spin states through the following assignments: 



£ T = Y + 

Sq Aq j 

£q Xq 1 

T _ - 

Xq 1 

vi = 



Xq ■ 



(A.12) 
(A.13) 
(A.14) 
(A.15) 



In the relativistic limit the neutrino momentum spinors 
become 



w(q,T) 
w(q,l) 







T ' 







(A.16) 
(A.17) 



while the antineutrino momentum spinors become 



«(q,T) 



(A.18) 
(A.19) 



As mentioned above, the fact that neutrino interactions 
involve only left-handed neutrino fields Pl f(y) implies 
that to 0(m 2 /E v ) only negative-helicity neutrinos and 
positive-helicity antineutrinos are produced. 
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Consider the density function defined in Eq. (|38|) ; see 
also Eq. Applying Eqs. (|A.2p and (|A.3p to neutrino 

fields with mass indices i,j, it can be expressed 



where 



Ni m (q,u) 



*C(q,u) = 



d 3 q d 3 u 

(2^r)3 (2^r)3 

- e'^-^N^q,!!) 



i(«j-*-«-»)N^»(q,u) 
(A.20) 



^^ ^ u ' (q,r,i) "' B(u ' M) 

X<4,»jOq,r,i>, (A.21) 

7kvi:^^ (q ' r '^ m(u ' s ' j) 



In these expressions the four-momentum (u^) = (u°, u) 
should not be confused with the momentum-space Dirac 
spinors u e (q, r, i). The components of the Pauli conju- 
gate spinors are u m = X)n( w *)"(7°)™ m an( ^ similarly for 
v m . 

To 0{ml/E v ) only a single 2x2 block of the 4 x 4 
spinor-space structure remains nonzero, and only one 
of the two spin states contributes to the expectation 
value. In particular the nonzero 2x2 blocks N^-- R (q, u) 
and N^(q, u) are those that would be projected out if 
(Nf™(q, u)) and (Nf™(q, u)) were sandwiched between 
the left- and right-projection matrices Pl and Pr. Em- 
ploying Eqs. fQ7)) and (|A~T8)) in Eqs. |A~2T|) and (TA~22|) 
results in the expressions 

N&*(q,u) - ^^(fl^/,,!,,), (A.23) 



N 



„ (q,u) - ^^Mft^ftu.tj) (A.24) 



L.R/ 



u = q, by virtue of the identity 



Sq Sq '/q '/q 



(A.25) 



with |q| = i? q in the relativistic limit. (This identity 
follows from an explicit expression for x~ that satisfies 
Eq. (lA~Til : 



cos 



(f) 



(A.26) 



where the polar angle 9 and azimuthal angle <f> give the 
direction of q.) Hence 



(A.27) 



for these non-zero 2x2 blocks. A nice form results when 



N^(q,q) = |^<6 q , T / q ,T,i) (A.28) 

in the relativistic limit with u = q. The trace of these 
2x2 blocks is 

Tr[N^(q,q)] = (< u a q ,u>, (A.29) 
Tr[Nf/(q,q)] = <^ T / q , Tj ). (A.30) 

That the leading factor becomes unity is a consequence 
of the tracelessness of the Pauli matrices cr. 
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